Non-Newtonian fluid flow through three-dimensional disordered porous media 
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We investigate the flow of various non-Newtonian fluids through three-dimensional disordered 
porous media by direct numerical simulation of momentum transport and continuity equations. 
Remarkably, our results for power-law (PL) fluids indicate that the flow, when quantified in terms 
of a properly modified permeability-like index and Reynolds number, can be successfully described 
by a single (universal) curve over a broad range of Reynolds conditions and power-law exponents. We 
also study the flow behavior of Bingham fluids described in terms of the Herschel-Bulkley model. 
In this case, our simulations reveal that the interplay of (i) the disordered geometry of the pore 
space, (ii) the fluid rheological properties, and (Hi) the inertial effects on the flow is responsible for 
a substantial enhancement of the macroscopic hydraulic conductance of the system at intermediate 
Reynolds conditions. This anomalous condition of "enhanced transport" represents a novel feature 
for flow in porous materials. 

PACS numbers: 47.56.+r, 64.60.ah, 47.50.-d, 47.11.-j 



The research on flow through porous media has great 
relevance for many problems of practical interest in sev- 
eral fields, including physics, medicine, biology, chemical 
and mechanical engineering and geology [l|, [SJ, HI ■ The 
disordered aspect of most natural and artificial porous 
materials is directly responsible for the presence of lo- 
cal flow heterogeneities that can dramatically affect the 
behavior, for example, of the transport of heat and mass 
through the system. Under this framework, the standard 
approach to investigate single-phase flow in porous media 
is to apply Darcy's law |l|, |2|, La] , which simply assumes 
that a global permeability k relates the average fluid ve- 
locity Uq in the field with the pressure drop Ap measured 
across the system, 
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where L is the length of the sample in flow direction and p 
is the viscosity of the fluid. As a macroscopic index, the 
permeability reflects the relation between the complex 
pore space morphology and fluid flow. 

In previous studies 0, IE 0, 0, S, S E3, detailed 
models of pore geometry have been used in combination 
with computational fluid dynamics simulations to pre- 
dict permeability coefficients and validate classical semi- 
empirical correlations for real porous materials. In princi- 
ple, the original concept of permeability as a global index 
for flow in porous media, however, is only applicable in 
the limit of Stokesian flow (linear). Strictly speaking, the 
validity of Darcy's law should be restricted to (i) New- 
tonian fluids and (ii) flows under viscous conditions, i.e., 
flows at very low Reynolds numbers, defined usually as 
Re = puodp/fi, where p is the density of the fluid and d p 
is the grain diameter. The departure from Darcy's law 
due to the contribution of inertial forces (convection) to 
the flow of Newtonian fluids has been the subject of sev- 
eral studies in the past 1(J U, 13] • In particular, it has 




FIG. 1: (Color online) Non-Newtonian (power-law) fluid flow 
through a typical realization of the Swiss-Cheese pore space 
(e = 0.7). The fluid is pushed from left to right. The solid 
lines with arrows correspond to trajectories of tracer particles 
released in the flow, while the contour plots give the velocity 
magnitude at different cross-sections of the porous medium. 
Their colors ranging from blue (dark) to red (light) correspond 
to low and high velocity magnitudes, respectively. 



been experimentally and numerically observed that the 
breakdown of condition (ii) can take place even under 
laminar flow conditions, i.e., before fully developed tur- 
bulence effects become relevant to momentum transport. 

In order to understand the physics of important prob- 
lems like blood flow through the kidney [13| or oil flow 
through porous rocks |l4j |. for example, one has to over- 
come the restriction (i) mentioned above by explicitly 
considering the nonlinear behavior of these fluids un- 
der shear, namely, their specific non-Newtonian proper- 
ties. Although these fluids have been known for a long 
time, technological applications which directly make use 
of their anomalous rheological behavior have come into 
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focus only recently. For instance, shear thinning solvents 
are present in dropless paints [15| . and shear thickening 
fluids are currently used as active dampers and compo- 
nents of enhanced body armors [l^]. While the physi- 
cal properties of Newtonian fluid flow through irregular 
media are theoretically well understood and have been 
confirmed by many experiments, non-Newtonian systems 
17 , HI, [Til lack a generalized description. In this Letter 
we investigate the flow of non-Newtonian fluids through 
three-dimensional porous media by direct numerical sim- 
ulation of momentum and continuity equations. To the 
best of our knowledge, this is the first time that non- 
linear effects coming from both rheological and inertial 
aspects of the fluid flow are considered simultaneously 
in the framework of a disordered three-dimensional pore 
space. 

The porous medium studied here is a three- 
dimensional realization of the Swiss-Cheese model [20j j . 
Spherical particles (solid obstacles) of diameter d p are 
sequentially and randomly placed in a box of length L in 
the x-direction and square cross-section of area A. Par- 
ticle overlap is allowed and the allocation process contin- 
ues up to the point in which a prescribed value for the 
porosity (void fraction) e is achieved. The mathematical 
formulation for the fluid mechanics in the interstitial pore 
space is based on the assumptions that we have a contin- 
uum and incompressible fluid flowing under steady-state 
and isothermal conditions. Thus, the momentum and 
mass conservation equations reduce to, 



pu-Vu = -Vp + VT 
V-u = , 
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where u and p are the velocity and pressure fields, re- 
spectively, and T is the so-called deviatoric stress tensor 
given by, 
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is the strain rate tensor. The 
variable fi is the dynamic viscosity, for which a constitu- 
tive relation must be provided in order to describe the 
specific non-Newtonian behavior of the fluid. Here we 
investigate the flow of two different types of rheologies, 
namely, the cross-power-law fluid and the Bingham fluid. 
The constitutive relation for a PL fluid can be written as, 



fll < /! < jU 2 , 



(5) 



where the constants [i\ and fj,2 are the lower and upper 
cutoffs, j = \/ \sijSij is the effective strain rate, K is the 



consistency index, and n is the power-law exponent. For 
n = 1 we recover the behavior of a Newtonian fluid. Flu- 
ids with n > 1 are shear-thickening, while shear-thinning 
behavior corresponds to n < 1. 

In the case of Bingham fluids, the rheology is com- 
monly approximated by the Herschel-Bulkley model 21 



22j which combines the effects of Bingham and power-law 
behavior for a fluid. For low strain rates, 7 < To/fj,o, the 
material acts as a very viscous fluid with viscosity fiQ. As 
the strain rate increases and the yield stress threshold tq 
is surpassed, the fluid behavior is described by, 
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where Kb is the consistency factor and n is the power- 
law index. Here we restrict our simulations to the case 
of Bingham fluids n = 1, i.e., the fluid is still Newtonian 
at large strain rates, with a viscosity [i = Kb- 

Non-slip boundary conditions are applied along the en- 
tire solid-fluid interface and end effects on the flow field, 
which become significant at high Reynolds numbers, are 
minimized by attaching ancillary zones at the inlet and 
outlet of the two opposite faces in the direction of the flow 
(i.e., ir-direction) . At the inlet, a constant inflow veloc- 
ity in the normal direction to the boundary is specified, 
whereas at the outlet we impose gradientless boundary 
condition. Finally, the four remaining faces are consid- 
ered to be solid walls. 

For a given realization of the porous medium and a 
given set of flow and constitutive parameters of the fluid, 
the numerical solution of the partial differential equations 
([2|) for the local velocity and pressure fields in the fluid 
phase of the void space, head and recovery zones is ob- 
tained by discretization using the control volume finite- 
difference technique 0. An unstructured grid with up 
to three million tetrahedral cells is adapted to the ge- 
ometry of the porous medium. For comparison, entirely 
consistent numerical solutions have also been calculated 
with a finite- volume scheme 24j. Finally, from the area- 



averaged pressures at the inlet and outlet positions, the 
overall pressure drop can be readily calculated. 

In Fig. [1] we show a three-dimensional plot of a typ- 
ical realization of the porous medium through which a 
power-law fluid flows. Clearly, the complex geometry of 
the pore space induces preferential channels on the flow 
whose localization and strength are significantly depen- 
dent on the rheological properties of the fluid as well as 
on the imposed inlet-outlet boundary conditions. For PL 
fluids, this intricate interplay between geometry and flow 
can nevertheless be macroscopically quantified in terms 
of an analogous to a permeability index, namely a hy- 
draulic conductivity, defined in terms of Darcy's law as 
kr> = K-yUoLJ Ap, where K\ is a reference viscosity taken 
as the consistency index for n = 1. As shown in the in- 
set of Fig. [2J the general behavior of kr> is qualitatively 
similar for different values of the exponent n. Moreover, 
it follows the characteristic trend of a simple Newtonian 
fluid (n — 1), namely that kr> remains essentially in- 
variant for low Re values up to a crossover point Re x 
where it starts to decrease due to the onset of non-linear 
convective effects on the flow 0, EH H^] . Quantitatively, 
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FIG. 2: (Color online) Flow of power-law fluids through three- 
dimensional porous media. The inset shows the variation of 
the ratio ko/k r with Reynolds number Re = puod p /Ki for 
different values of the power-law exponent n and e = 0.5. The 
resulting data collapse presented in the main plot confirms the 
adequacy of our rescaling procedure in terms of the modified 
permeability index k n /ko and the modified Reynolds number 
Re n (see text). 



however, we observe that both the upper limit for ku and 
Re x are strongly dependent on n. 

Darcy's law has been gen eralized to power-law fluids in 
previous studies 2(| 27 , 28| . Here we define an hydraulic 
conductivity as, 
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As shown in Fig. [3J this generalized index when calcu- 
lated at low Reynolds numbers, namely k , can be consis- 
tently correlated with intrinsic properties of the fluid and 
porous medium by means of the following semi-empirical 
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expression 
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where d e is the only fitting parameter corresponding to 
an average effective pore diameter, namely, the average 
pore size (in units of d p ) of the system calculated as if it 
was a packed bed consisting of identical spheres 2(| 2i| . 
The parameter k r corresponds to the value of fcrj cal- 
culated for a Newtonian fluid (n = 1) under very low 
Reynolds conditions, i.e., the porous medium permeabil- 
ity according to Darcy's law. 

In order to substantiate the non-Newtonian aspect of 
the fluid, it is also necessary to redefine the Reynolds 
number as [29(, 
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FIG. 3: (Color online) Dependence of the hydraulic conduc- 
tivity at very low Reynolds numbers ko on the power-law ex- 
ponent n for two different values of porosity e. The solid lines 
are the least-squares fits to the simulation data using Eq. © 
with d e /dp — 0.35 and 1.58, for e = 0.5 and 0.7, respectively. 



where the term (1 — e)/e 3 has been adapted from the 
classical Kozeny-Carman equation [l|. It is worth men- 
tioning that Eq. ([9]) breaks down close to the critical 
percolation porosity [30( . 

In the main plot of Fig. [2] we show that all data sets 
of k n /ko against Re„, with fco obtained from Eq. ©, 
collapse onto a single curve for the entire range of (mod- 
ified) Reynolds numbers, independent of the numerical 
values of e and n. Despite the details of the model porous 
medium geometry employed here as well as the complex- 
ity of the fluid rheology, this remarkable invariance of 
behavior suggests that the resulting flow properties of 
the system remain in the same universality class of New- 
tonian fluid flow in disordered porous media. 

Next we present results for flow through three- 
dimensional porous media of Bingham fluids with rheol- 
ogy given approximately by the Herschcl-Bulkley model 
Eq. ([6]). The proper way to quantify inertial and vis- 
cous forces in this case is to define the Reynolds number 
as Res = pu d p /KB- In Fig. 0] we show that the lin- 
ear hydraulic conductivity, defined as ks = KbUqL/ Ap, 
remains constant up to a certain crossover that is pro- 
portional to the threshold To, Re x ~ tq. Below this 
crossover, since the fluid has Newtonian behavior with 
high viscosity /Uo everywhere in the pore space, the flow 
can be macroscopically described by Darcy's law. Above 
this crossover, the presence of low and high strain rates 
zones in the flow leads to a nonuniform spatial distri- 
bution of fluid viscosity, therefore increasing the over- 
all permeability index ^ of the system. This behavior 
persists up to the point in which inertial forces become 
relevant. While the specific fluid rheology investigated 
here tends to enhance the flow at high Res, the effect 
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lems, e.g., chemical reactors, chromatographic columns 
and switches for flow. Finally, at sufficiently large values 
of Res , the viscosity of the fluid is uniform and therefore 
the local permeabilities become all the same, regardless 
of the value of Res and tq. In this situation, all curves 
of ks collapse. 

Summarizing, in spite of the non-linear nature of the 
fluid rheology and the complex geometry of the inter- 
stitial pore volume, in the case of power-law fluids, we 
have shown the remarkable fact that the flow behavior 
can still be quantified in terms of an universal curve ex- 
tending over a broad range of Reynolds conditions and 
power-law exponents. Our results for Bingham fluids are 
even more striking. There the pore space geometry, fluid 
rheology and inertia can combine to generate a partic- 
ular condition of "enhanced transport" which should be 
found in experiments. 

We thank CNPq, CAPES, FUNCAP, FINEP, Petro- 
bras, the National Institute of Science and Technology 
for Complex Systems , and the Swiss National Science 
Foundation (SNF) under Grante No. 116052 for finan- 
cial support. 




FIG. 4: (Color online) Flow of Bingham fluids (Herschel- 
Bulkley model) through three-dimensional porous media. The 
plot shows the variation of the ratio ks / k r with Reynolds 
number Res for different values of the parameter to, as de- 
fined in Eq. (B). Here k r corresponds to the lower limit of fcs 
at very low Reynolds conditions. The presence of maxima in 
all cases is a distinctive result of the competition between rhe- 
ology and convective non-linearities. The contour plots in (a), 
(b) and (c) show the spatial variation of the magnitude of the 
local ratio |w| / |Vp| calculated on the cross-section through 
the middle of the porous medium parallel to the flow, for 
r = 0.1 and Res = 3.5 x 10~ 2 , 1.7 and 35, respectively. Their 
colors ranging from blue (dark) to red (light) correspond to 
low and high values of \u\ / |Vp|, respectively. 



of inertia is to reduce the permeability index under the 
same conditions 0, E3] • As a result of this competition, 
a maximum hydraulic conductivity can be observed at 
an intermediate value of Res that is also dependent on 
the threshold To. As shown in Fig. 2J this effect is bet- 
ter illustrated when we observe contour plots of the local 
ratio \u\ / |Vp| calculated at the middle cross-section of 
the porous medium. To the best of our knowledge, this 
condition of "enhanced flow" through disordered porous 
media represents a novel regime of momentum transport 
that could have potential applications in practical prob- 
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